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1 Introduction 



At high temperatures the primeval cosmic plasma is very close to equilibrium. This 
is because the reaction rate, T = an, typically is much higher than the expansion 
rate, H ~ 1/1 Here a is the characteristic cross-section and n ~ T 3 is the particle 
number density. One can see from the covariant Boltzmann kinetic equation: 

(d t -Hpd p )f = I coll (1) 

that equilibrium for massless particles is not distorted by the expansion of the uni- 
verse. Indeed the equilibrium distribution: 

f eq =[ex P (E- t i(t))/T(t)±l]- 1 (2) 

which annihilates the collision integral I co u, satisfies this equation if T/T = —H 
and fi(t) ~ T(t). However if the particle mass m is nonzero, the l.h.s. of eq. (0) 
cannot vanish for all possible values of the momentum p with any choice of the 
functions T(t) and The deviation from equilibrium can roughly be estimated as 
5f/f = Hm 2 /(TTE). 

Since photons and electrons are very strongly coupled, the deviations from equi- 
librium in the electromagnetic cosmic plasma are tiny. This is the reason for a perfect 
Planckian shape in the observed spectrum of the cosmic microwave radiation, where 
the distortion is known to be smaller than 10 -4 . One would expect the same to be 
true for massless cosmic neutrinos. However this is not so and the spectral correc- 
tions, especially in the high energy tail, can be of several per cent. The physics of 
the phenomenon is quite simple. At high temperatures neutrinos and electrons are 
sufficiently strongly coupled to each other, they have the equilibrium distributions 
and their temperatures are equal. At smaller temperatures, below 2-3 MeV, the neu- 
trinos decouple and below these temperatures the primeval plasma consists of two 
(or better to say four) almost decoupled components, the electromagnetic and three 
neutrino ones. In the approximation of instantaneous decoupling the neutrino spec- 
tra maintain their equilibrium shape with the temperature decreasing as the inverse 
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scale factor, T v ~ l/a(t). The temperature of electrons and photons decreases slower 
because e + e~-annihilation heats up the electromagnetic part of the plasma when the 
temperature drops below m e . Ultimately the ratio of the temperatures reaches the 
well known value T 7 /T^ = (11/4) 1 / 3 = 1.401, see e.g. ref. Jl|. This result was ob- 
tained under the assumption of entropy conservation in the electromagnetic sector, 
which is not necessarily true in the non-equilibrium case, where the limiting value of 
T 7 /T^ becomes slightly different. The interactions between the electromagnetic and 
neutrino components of the plasma, though weak, are not completely vanishing and 
therefore the residual annihilation e + e~ — > vv and, to a lesser extend, elastic scatter- 
ing ev — > ev slightly heat up the neutrinos and distort their spectra. The distortions 
are larger at higher energies because the weak interactions are getting stronger with 
rising energy. 

The overall neutrino heating under assumption of equilibrium spectra has been 
considered earlier in refs. ||, || ||. In this approximation the problem is very much 
simplified but the method itself is far from being precise from the very beginning. 
A more complete approach with the distortion of the spectra taken into account has 
been discussed in several papers |5|, ||, [jj in different approximations and with some- 
what different results. In what follows we present an accurate numerical treatment 
of the problem with a better accuracy than that of the preceding paper where 
the numerical solutions of the exact kinetic equations were also obtained. Our bet- 
ter accuracy is achieved by a convenient reduction of the collision integral to two 
dimensions and a powerful integration method developed in ref. ||. 

2 Basic Equations 

Instead of time and momenta we choose the following dimensionless variables: 

x = ma(t), i/j = Pjd(t) (3) 
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where m is an arbitrary parameter with dimension of mass, which we took as m = 
1 MeV and the scale factor a(t) is normalized so that a(t) = 1/T„ = 1/T 7 at high 
temperatures or at early times. In terms of these variables the kinetic equation (|I|) 
can be rewritten as: 

Hxd x f(x,y 1 )=I coU . (4) 

The collision integral I co u is dominated by two-body reactions 1 + 2^3 + 4 and is 
given by the expression: 

1 ^ f d 3 p2 d 3 p% d 3 p4 



l coll 



E 



2E l ^ J 2£ 2 (2tt) 3 2E 3 (2tt) 3 2£ 4 (2tt) 3 
(2n)W\ Pl + P2 -p 3 - p^Fih, / 2) / 3 , U)S | A\j 2 ^ u (5) 

where F = f 3 h{l ~ /OC 1 ~ /a) - fM^ - fa) (I ~ A), \ A \ 2 is the weak interaction 
amplitude squared summed over spins of all particles except the first one, and S is 
the symmetrization factor which includes 1/2! for each pair of identical particles in 
initial and final states and the factor 2 if there are 2 identical particles in the initial 
state; the summation is done over all possible sets of leptons 2, 3, and 4; The relevant 
reactions are presented in Table 1 for the case when the first particle is i/ e , and in 
Table 2 for or v T . Our results in these tables agree with those of ref. 0] except 
for the reactions v a v a — > v a v a where we took twice larger contribution because of 
identical particles in the initial state. 

We assume that the distribution functions for and v T are equal while that for 
v e is different because of different strength of charged current interactions. Similarly 
we assume that the lepton asymmetry is negligible so that f v = fp. Therefore there 
are two unknown functions of x and y: f Ve and f v = f Vr . Each of them satisfies 
the kinetic equation (^) with the appropriate collision integral. We also assume that 
the distributions of photons and electrons are the exact equilibrium ones (0) with 
the unknown temperature T 7 (x) and zero chemical potential, = 0. The third 
necessary equation is the covariant energy conservation: 

sfr+n (6) 



where p is the total energy density: 



vr 2 T 7 4 2 r dqq 2 ^Jq 2 + m 
P ~ ~1E~ + ^ 2 J exp (E/T y ) + 

and P is the pressure: 



r aqq \ q + m t l r 2 r 



45 n 2 J 3yV + m 2 [exp(£/T 7 ) + l] 3tt 2 7 3tt 2 7 

The Hubble parameter, H = a/ a, which enters the kinetic equation (^) is ex- 
pressed through p in the usual way, 3H 2 m 2 Pl = 8irp, ignoring the curvature term and 
the cosmo logical constant. 

The collision integral in eq. (pj) can be reduced from nine to two dimensions as 
described in Appendix A. For the sake of brevity let us introduce some notations: 
f a ( Pj ) = /0'>, d x = D u 4(3,4) = D 2 (S,A)/E 3 E 4 , and 4 = D 3 /E X E 2 E 3 E A . The 
functions Dj are defined in Appendix A. We do not indicate the arguments for the 
functions Di and D 3 because they are symmetric in all 4 arguments. In terms of these 
functions we can write the coupled kinetic equations for f Ue and f v in the following 
way: 

G 2 r 

Hxd x fl 1 } = J dp 2 p 2 dp z p z dp 4 p 4 b{E x + E 2 - E 3 - E 4 ) 

{F[f£\ f£\ /i 3) , /«] [64 - 4d 2 (l, 4) - 44(2, 3) + 24(1, 2) + 24(3,4) + 64] 
+Fif£\ 4 2) , f£\ /«] [44 + 2rf 2 (l, 2) + 2rf 2 (3, 4) - 24(1, 4) - 24(2, 3) + 44] 

+F[f£\fg\ fl% /«][2dx - 24(2, 3) - 24(1, 4) + 24] 
+F[f£\ / e (2) , /i 3) , / e (4) ] W L + si)(24 - 4(2, 3) - 4(1, 4) + 4(3, 4) 
+4(1, 2) + 24) - % L ^m 2 (4 - 4(1, 3))/£ 2 £ 4 ] 
+F[fit\ /i 2) , / e (3) , / e (4) ] [4^1(4 - 4(2, 3) - 4(1, 4) + 4) + 
Ag 2 R (4 - 4(2, 4) - 4(1, 3) + 4) + 4^m 2 (4 + 4(1, 2))/E 3 eA (9) 



and: 



G 2 r 

Hxd x f^ = —j— J dp 2 p 2 dp 3 p 3 dp i p 4 5{E 1 + E 2 - E 3 - E 4 ) 



2n 3 pi 



{F[/g\ f®,f$>, /«] [9(1, - 6rf 2 (l, 4) - 6^(2, 3) + 34(1, 2) + 34(3, 4) + 94] 
/«, /W , /«] [24 + 4(1, 2) + 4(3, 4) - 4(1, 4) - 4(2, 3) + 24] 

/«,/£), /£>] [4 - 4(2, 3) - 4(1, 4) + 4] 

+^ ) , /i 2) , /i?, /i 4) ] [4(^1 + aD(2d! - 4(2, 3) - 4(1, 4) + 4(3, 4) 

+4(1, 2) + 24) - shgnmlid! - 4(1, 3))/£ 3 £ 4 ] 

+F\f$, / e (3) , / e (4) ] [4^1(4 - 4(2, 3) - 4(1, 4) + 4) + 
4^(4 - 4(2, 4) - 4(1, 3) + 4) + 4^77^(4 + 4(1, 2))/E 3 E 4 }} , (10) 

where g L = 1/2 + sin 2 9 Wl g R = sin 2 9 W and g L = —1/2 + sin 2 9 Wl f e = [exp(E/T J ) + 
and the temperature T 7 is determined from eq. (|6|). Thus we have the com- 
plete system of three equations for three unknown functions T 7 (x), f Ue (x,y), and 
fva i x -> v) which we solve numerically. The solution has been found in two different 
but equivalent ways. First, the system has been solved directly, as it is, for the full 
distribution functions f v . (x, y) and, second, for the small deviations 5j from equilib- 
rium f Uj (x,y) = f(f(y)(l + 5 3 (x,y)), where fM = [exp{E/T v ) + 1]' 1 with T v = l/o. 
In both cases the numerical solution has been done exactly, not perturbatively. So 
with infinitely good numerical precision the results must be the same. However since 
the precision is finite, different ways may give different results and their consistency is 
a good check of the accuracy of the calculations. It is convenient to introduce S(x,y) 
because the dominant terms in the collision integrals, containing neutrinos only, can- 
cel and subdominant terms are proportional to S. In the parts of the collision integrals 
which contain electrons, there are also contributions proportional to the difference in 
temperatures (T 7 — T v ). We have found that the results of these two methods, i.e. 
solving either for f u or 5, are very close (see the next Section and Table 3). 

To check possible algebraic errors we have made the following procedure. First, we 
have solved numerically the three relevant equations in the Boltzmann approximation, 
simply changing the function F which is defined after eq. (P) to F — > f%f\ — /1/2 
with the same coefficients and d-functions in the collision integral as presented above. 
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Second, we repeated the numerical calculation but this time we changed the terms 
describing the direct reactions (the ones that enter I co u with negative signs) with 
(presumably) identical terms where the integration over the phase space of the final 
particles had been made explicitly. That is, the positive terms were described through 
d-functions as above, while negative ones had been calculated in a different and 
independent way. The expressions for the negative terms are presented in Appendix B. 
The results of these two numerical calculations are in perfect agreement. This excludes 
possible accidental errors in calculations of D-functions (see Appendix A) as well as 
in the typing of the program for the numerical integration. 

3 Description of Calculations and Results 
3.1 Initial conditions 

We have numerically calculated the evolution of the distribution functions satisfying 
the system of eqs. (|||9],[10]) in terms of the dimensionless time x = m^a in the interval 
Xi n < x < Xf, where a is the scale factor. With our normalization the "time" x = 1 
corresponds to the neutrino temperature T u = itlq = lMeV^\ 

For large values of x the collision integrals in the r.h.s. of eqs. (|9|,|T0D are sup- 
pressed by the factor 1/x 2 . We have found that at x ~ 50 all the functions f Ue , 
/„ m(t) , and T 7 reach their asymptotic values, so to be certain that we are in the 
asymptotic limit we have chosen the final time Xf = 60. For the initial time Xi n 
we have chosen three different values Xi n = 0.1, 0.2, and 0.5 which correspond to 
T u = lOMeV, 5MeV, and 2MeV respectively. In all the cases we have assumed that 
the neutrinos are in thermal equilibrium for x < X{ n with the distribution functions 
/if) = l/[exp(y) + l]. 

For the dimensionless momentum y = pa we took a 100 and a 200 point grid, 

equally spaced in the region < y < 20. In contrast to ref. we have included the 

5 Strictly speaking for non-equilibrium neutrinos we cannot introduce temperature T v . But since 
non-equilibrium corrections to the neutrino spectrum are small (see below) we sometimes use the 
notation T v instead of 1/a, because in equilibrium T v a = 1. 
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point y = in which the collision integrals in eqs. (|9l, |l0|) have different analytical 
expressions due to the cancellation of the factor l/p\ in front of the collision integrals 
by the corresponding factors in the functions D{. This allows us to compare the 
collision integral at y = with the integral in a nearby point < y <C 1 in order to 
check that the numerical errors are small in the region of small momenta y < 1. 

3.2 Numerical results 

There are three phenomena which play essential roles in the evolution of the neu- 
trino distribution functions. The first of them is the temperature difference between 
photons and e on one side and neutrinos on the other, which arises due to the heat- 
ing of the electromagnetic plasma by e + e~-annihilation. Through the interactions of 
neutrinos with electrons this temperature difference leads to non-equilibrium distor- 
tions of the neutrino spectra. The temperature difference is essential in the interval 
1 < x < 30. The second effect is the freezing of the neutrino interactions because 
the collision integrals in the r.h.s. of eqs. fl9[,|10D drop as 1/x 2 . At small i<1 the 
collisions are very efficient but at x > 1 they are strongly suppressed. The third 
important phenomenon is the elastic ////-scattering which smoothes down the non- 
equilibrium corrections to the neutrino spectrum. It is especially important at small 
x < 1. 

Accordingly, the evolution of the distribution functions proceeds as follows. At 
x < 0.2 the temperature difference is negligibly small while ////-scattering is strong. 
Therefore the neutrino distribution functions are practically in equilibrium. In the 
interval 0.2 < x < 4 all three effects are essential and during this period the neutrino 
spectra are distorted and simultaneously the photon temperature slightly drops down 
because of the energy transfer from to neutrinos, as compared to the photon 
temperature in the approximation where the neutrino decoupling is instantaneous. 
When x > 4 the collision integrals in the r.h.s. of eqs. (|9|JT0|) are small and the 
temperature difference practically does not affect the neutrino distribution functions, 
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even though the temperature difference is the largest in this time interval. 

In Figs. 1-5 we present the results of the numerical solution of the system of 
eqs. (^yiDp in the interval 0.1 < x < 60 (approximately 4000 points in x) both for 
Fermi-Dirac (FD) and Maxwell-Boltzmann (MB) statistics. In these runs we used a 
100 point grid in the momentum interval < y < 20. 

The ratio T^jT v as a function of x is presented in Fig. 1 for FD and MB statistics. 
If we neglect the exchange of energy between and neutrinos in eq. (|6|) (which 
is equivalent to entropy conservation) the asymptotic values are T^jT v = 1.4010 for 
FD statistics and T 1 jT v = 1.4423 for MB statistics. When the energy exchange is 
taken into account (energy conservation) these ratios have slightly different values 
T^/Ty = 1.3991 (or 1.3994, using the full distribution function) for FD statistics and 
Tj/% = 1.4404 for MB statistics. 

In Fig. 2 the correction to the total neutrino energy density bp v j p eq for the case 
of FD statistics is plotted. Here p eq = (77r 2 /240)(m/x) 4 is the unperturbed energy 
density of neutrinos and: 

OPu = g - Peg , (11) 

with p Ue and p v found from our numerical solution. The upper curve in Fig. 2 
corresponds to the entropy conservation case, while the lower one corresponds to 
the energy conservation case. The effect is smaller in the second case because the 
"driving force", which is proportional to the temperature difference, is lower due to 
the excessive photon cooling. We can see that while the difference in the asymptotic 
values of T 1 jT v for these two cases is as small as 0.14% the difference in 5p v / p eq is as 
large as 20%. This can be explained in the following way. Compare Fig. 1 and Fig. 
2 for the FD case. We see that already at x = 4 the value of the function 5p u / p eq 
is close to its asymptotic value. At the same moment, x — 4, the ratio T^/T v has 
a value of 1.08 which is much smaller than the asymptotic one, 1.401. This means 
that the dominant contribution to the distortion of the neutrino spectra comes from 
the period when the temperature difference AT = T 7 — T v was rather small. The 
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smallness of AT was compensated by a more efficient energy exchange between v and 

at smaller x. In this range of x, a relatively small variations in T 7 /T, connected 
to different approximations made in the calculations, should be compared with the 
small difference AT. Thus a small variation in T 7 /T, is relatively enhanced, and we 
see how these small variations can produce the strong 20% effect. 

In Fig. 3 the relative deviations of the neutrino spectra from the equilibrium one, 
5j(x,y) = (f u . — f eq )/feq, are presented as functions of x for several values of the 
momentum y = 3,5, 7. In Fig. 3a we present the results for electronic neutrinos 
and in Fig. 3b for muonic (tau) ones. For large momenta y the deviation from the 
equilibrium is larger and the maximum asymptotic value is reached later (at larger 
x). This is a result of stronger interactions of more energetic neutrinos. 

In Fig. 4 we compare the deviations from the equilibrium distributions, 5 Ue and 
5j/ . , for FD and MB statistics. We plot 5i for the same value of the momentum 
y = 5 as functions of x. We see that the results for the case of Boltzmann statistics 
are larger than those for the Fermi one by approximately 25%. For both FD and MB 
statistics the spectral distortion for u e is more than twice larger than that for or 
v T . This is due to a stronger coupling of v e to e ± . 

In Fig. 5 we plot the asymptotic values of the corrections to the neutrino distribu- 
tions 5j = (f Uj — ft q )l ft q as functions of the dimensionless momentum y. The dashed 
lines a and c correspond to Maxwell-Boltzmann statistics and the solid lines b and 
d correspond to Fermi-Dirac statistics. The upper curves a and b are for electronic 
neutrinos and the lower curves c and d are for muonic (tau) neutrinos. All the curves 
can be well approximated by a second order polynomial in y like 5 = Ay(y — B) in 
agreement with ref. ||. 

3.3 Checking the results 

First, let us discuss the choice of the initial time X{ n . We made the runs for the system 
of eqs. (|5|,P|,|lTi|) with three different values x in = 0.1,0.2 and 0.5. We found that the 
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results of the runs with Xi n = 0.1 and X{ n = 0.2 are the same with the accuracy of 
10~ 5 in the distribution function for all values of x (see Table 3). This means that for 
x < 0.2 we can neglect the non-equilibrium corrections to the neutrino distribution 
functions. For larger x this is not so and e.g. the spectral distortion calculated with 
Xi n — 0.5 is smaller by approximately 15% than the one calculated with Xj„ = 0.2. 

In order to check the errors connected with a finite number of points in the mo- 
mentum interval < y < 20 we took a 100 and a 200 point grids. We found that 
the difference in the results is small, which means that a 100 points grid is a good 
approximation for the numerical solution of our system of kinetic equations. 

The errors connected with a finite number of points in time x are much more 
important. We control these errors in the following way. First, we run the program 
with some fixed number of points in time x, distributed in the time interval Xi n < x < 
Xf in such a way that the distribution functions do not change significantly at any 
momentum point y during one time iteration dx. Then we run the program for the 
entropy conservation law (i.e. with equilibrium neutrinos) with the same values of 
time Xi as in the first run. Then we compare the asymptotic values of the temperature 
ratios with the sample values which are T^/T v = (11/4) 1//3 = 1.4010 for FD statistics 
and T 7 /Tj, = 3 1/3 = 1.44225 for MB statistics. We found that when we run with 400 
time points, the numerical error in these temperature ratios are of the order ~ 0.001, 
and this value is only one order of magnitude smaller than the effect itself, which is 
not good enough. However, in the case of 4000 time points this error is as small as 
~ 10~ 4 . In our calculations we used a 4000 point grid in time. 

We also checked the errors related to the finite number of points in time in a 
different way. Instead of the simple time evolution we used the Bulirsch-Stoer method, 
described in the book |J. We found, that this method gives the same results as the 
simple time evolution with 4000 time points (see the first two rows in Table 3). 

Besides reducing the numerical errors we have used two different analytical ap- 
proaches in order to check the results. In the case of FD statistics we solve the system 
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of eqs. (|BtP,|lO|) both for the full neutrino distribution functions f Uj and for the de- 
viations from equilibrium Sj = (f Vj — f„ q )/ ' fl q . In the last case the contributions 
to the collision integrals for all the processes vanish for vanishing S, except for the 
interactions of neutrinos with electrons, where the "driving force" term, proportional 
to the temperature difference between v and e , gives a nonzero contribution. Com- 
paring the results in these two cases (with the same values of all parameters) we have 
found that the non-equilibrium corrections to the neutrino spectra are systematically 
smaller by about 10% when the equations are solved for the full distribution functions 
f v . . The comparison of the results of the different ways of calculation are summarized 
in Table 3. 

3.4 Helium abundance 

We have modified the standard nucleosynthesis code (ref. ]10|) in the following way. 
The neutrino distribution functions have been parametrized as: 

fu^U-^ + SiipM), (12) 

where 5, were fitted as functions of momentum and temperature in accordance with 
our numerical solutions. 

First, the neutrino energy densities p Ui have been calculated with the corrected 
spectra. This changes the universe cooling rate and has the major impact on the 
helium abundance. 

Second, the 6 weak interaction rates for (n <-» p)-reactions have been modified with 
the account for the spectral distortion of the electronic neutrinos. The 3 reactions 
with v e or v e in the initial state, n v e — > p e~, p v e — > n e + and p v e e~ — > n, are 
more important and the neutrino distribution functions have been corrected according 
to eq. (0). However, for the 3 inverse reactions with neutrinos in the final state, 
n e + — > p i>, n — >• p e~ v and p e~ — > n v, the spectral modification results only in a 
small correction to the Pauli suppression factor, (1 — /„) — > (1 — f u ■ (1 + 5)), and the 
effect of the corrections to these last 3 reactions is negligible. 
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Finally the law of the evolution of the photon temperature is slightly changed, 
which in turn also affects the reaction rates. The correct photon temperature is 
calculated by including the neutrino energy densities into the equation governing the 
variation of the photon temperature in the course of the universe expansion: 



dT 7 p EM + p EM + 4/3 p v 



(13) 



d InV dpEM/dT^ + dp v j dT^ ' 

where V = a 3 . If neutrinos had the equilibrium distribution their energy density 
would satisfy the covariant conservation law (^) separately so that their contribution 
into eq. (|T3|) would cancel out. However this is not true if the energy exchange 
between z/s and e ± is taken into account. 

The final change in helium-4 abundance has been found to be AY = 1.4 • 10~ 4 . 
Similarly we find relative changes in the 3 He and 7 Li abundances of the order 5 • 10~ 4 
and —6 • 10~ 4 respectively. 

4 Discussion 

Our results agree reasonably well with the previous numerical calculations of ref. [[7]. 
For example we obtain 5p u Jp Ue = 0.94% (0.83%) and 6p v Jp Vfl = 0.40% (0.33%) (see 
Table 3), to be compared with 0.83% and 0.41% respectively obtained in ref. 0. 
For the shape of the spectral distortion the difference is larger, especially for low 
y. To characterize the difference let us introduce the effective neutrino temperature 
as T e jf = E/hn[(l — /)//], which we compare to the unperturbed temperature T 
calculated in the limit of entropy conservation i.e. in neglect of the neutrino heating. 
The difference: 

1kiz2± = l^l 5{x , y) (14) 

T y 

goes to a small negative value for small y in our case and to a considerably larger 
positive one according to ref. . We ascribe this difference to a better accuracy of our 
calculations which is most profound at low y (in particular we took a 100-200 point 
grid in comparison to 35 in ref. 0). The physical reason for this negative result 
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is the elastic z/e ± -scattering. Since electrons and positrons have a slightly higher 
temperature than the neutrinos, the scattering will depopulate the low momentum 
region pushing neutrinos to higher momenta. Our agreement with ref. |7| is much 
better for the physically more important region of high y. The distribution function 
decreases exponentially, whereas the correction S increases quadratically. Moreover 
the weak interaction rates are typically proportional to the energy squared, so the 
dominant contribution of 8 to neutron-proton reactions lays near y = 5. For large y 
our results are systematically larger than that of ref. |7| by approximately 15%. This 
difference cannot be ascribed to the difference by a factor two in the collision integrals 
for the scattering of identical neutrinos because that leads to the opposite sign of the 
effect and is rather small, giving corrections of a few per cent. We believe that the 
difference between us and ref. |7j arises because of different numerical accuracies. 

Another source of the difference in the results is the different treatment of the 
cooling of the electromagnetic plasma. The simplest approach is to neglect the loss 
of photon energy which goes into neutrinos. In this case the entropy of photons and 
e + e~-pairs is conserved and the function T 7 (x) can be determined from the relation 
x 3 {pem + Pem) /T~f = const. In this approximation the cooling of the electromagnetic 
plasma is induced only by the expansion of the universe. The energy transmitted 
to neutrinos results in an excessive cooling which we have calculated exactly from 
the more complicated equation (|6|), which does not necessarily imply conservation of 
entropy in the non-equilibrium case. The effect is found to be significant. Inclusion 
of the additional cooling of photons amplifies the neutrino spectral distortion by 
approximately 20% (see fig. 2). This can be explained as follows. The distortion of 
the spectrum is determined by the temperature difference AT/T U = (T y — T v )/T v and 
by the strength of neutrino-electron interactions. The former rises with the expansion 
while the latter dies down. At large x the temperature difference AT is much larger 
than the correction to the photon temperature but at smaller x = 0(1) they do 
not differ so much. The dominant contribution to the distortion of the neutrino 
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spectrum comes just from the period when x is near unity and AT is relatively small, 
so that the small correction to T 7 is relatively enhanced. This effect was neglected 
in the previous papers [||, [7| so their results are somewhat underestimated; it was 
approximately calculated in ref. || under the assumption that the neutrinos have 
equilibrium spectra. However the spectral modification gives a contribution into 
additional photon cooling of the same order of magnitude and should be taken into 
account. 

In the earlier papers ]5], |5| the effect was considered in the Boltzmann approxi- 
mation which simplifies the calculations very much. Another simplifying assumption, 
previously used, is the neglect of the electron mass in the collision integrals for ve- 
scattering and for annihilation vv — ■> e + e~. We have also made numerical calculations 
in these approximations and found that both of these give rise to a larger spectral 
distortion in agreement with ref. |7J]. In ref. |J the effect was calculated numerically 
while in ref. || an approximate analytical expression was derived. However in ref. 
H the influence of the back-reactions which smoothes down the spectral distortion 
was underestimated due to a numerical error in the integral. With the correction of 
this error the effect should be twice smaller (in the approximations of that paper). 
Our numerical results in the limit of Boltzmann statistics and for m e = are in a 
reasonable agreement with the corrected estimate of that paper: 

S Va ^3x lCr 4 y(lh//4-3). (15) 
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A Reduction of the collision integral 

In this Appendix we will make analytically as many integrations in the collision inte- 
gral as possible. In the case we are interested in, i.e. four-particle interaction vertex 
and the isotropic one-particle distribution functions, seven out of nine integrations can 
be done and only two integrals upon the momentum remain for numerical treatment. 
We start from the general form of the collision integral in the space-homogeneous 
case: 

'coll = (2-)V(E^)|M /l (|p|)| 2 F( / ) , (16) 



where the particle energy is E i = Jmj + p? and F(f) is defined as: 

F(f) = [i - mi - wj', - [i - /{][i - mf 2 , (I?) 

where the distribution functions f) only depend upon the moduli of the particle 
momentum and time fj(\p\,t). We make use of the identity: 

5 3 (Ep*) = / ei(A ' pl+P2 " p ' 1 ^ ) (0 ( 18 ) 

and explicitly separate out the angle integrations: 

d 3 pi = d<pid cos Oipfdpi = p^dpidVLi . (19) 
The collision integral in eq. ([H]) takes the form: 

'cC = ^S^^-^^^^^^ ■ <20) 
Here we have defined: 



L l>2 



J e-^dn p{ J e-«^dn p , 2 \M fi (\p\)\ 2 . (21) 

In the four-fermion approximation all the possible squared matrix elements consist 
of two kinds of terms: 

K l {q 1 ^){q^ql) = K X {E X E 2 - q x q 2 )(E 3 E 4 - qm) (22) 
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and: 



K 2 m 2 {q 3 <fi) = K 2 m 2 (E 3 E 4: - q 3 q 4 ) 



(23) 



where every qi corresponds to some of Pi. In order to perform the angle integrals we 
define the angles between <fj through angles from eq. (]2~]J) as: 



cos 0i2 = sin 6\ sin 6 2 cos(<^i — (p 2 ) + cos Q\ cos 6 2 ■ 



(24) 



The integration in the angles (pi or (p 2 over the total period and the structure of the 



matrix elements in eqs. ( p2| , p3|) result in the vanishing of the first term in eq. (|24]). 
Then all angle integrals in eq. (^TJ) can be trivially taken and we come to three kinds 
of integrals^: 

4 f°° dX 

Di = - — sin(Agi) sin(Ag 2 ) sin(Ag 3 ) sin(Ag 4 ) , (25) 

7T JO X Z 



D 2 (3, 4) = I ^ sin(Agi) sin(Ag 2 ) 



7T 



A 2 



cos(A<? 3 



sin(Ag 3 
Ag 3 



cos(Ag4) — 



sin(Ag4 



Xq 4 



(26) 



and: 



4gig 2 g3?4 r°° d\ 

n Jo A 2 

sin(Ag 3 



cos(Agi 



sin(Agi 
Agi 



cos(Ag 2 ) - 



sin(Ag 2 
Ag 2 



cos(Ag 3 ) - 



Ag 3 



cos(Ag4) — 



sin(Ag 4 
Ag 4 



(27) 



For the matrix element squared given by eq. ([22|) the function D in eq. fl21"|) can be 
written in the form: 

D = K 1 (E 1 E 2 E 3 E i D 1 + D 3 ) + K x [E X E 2 D 2 (Z,A) + E 3 E A D 2 (1, 2)] , (28) 



while for the matrix element squared given by eq. (123) we get: 



D = K 2 E X E 2 \E 3 E^D X + D 2 (3, 4)] . 



(29) 



6 In the case that one of or (74 in the second integral D2 corresponds to the incoming particle 
and the other to the outgoing one, D 2 changes sign. 
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The calculation of the integrals (|25| - [27| ) is straightforward. The functions D\ and 
D 3 are symmetric with respect to permutations of any variables and D 2 is symmetric 
under permutations 1 «-> 2 and 3 <-> 4. In what follows we present the expression for 
D2(3, 4). Without loss of generality we can assume that q± > q 2 and 93 > q^. Then we 
will have four different cases, which depend upon relations between the momenta.[] 

a) qi + q 2 > <?3 + <?4 and q 1 + q 4 > q 2 + q 3 , 

Di = ^{q 2 + q-i + qi-qi) (30) 
D 2 = ^(( qi - q 2 f + 2(q 3 3 + q 3 ) - 3( 9l - q 2 )(q 2 3 + ql)) , (31) 

^3 = ^ (ql ~ 5qfql + 5q 2 iq 3 2 - q 5 2 

- 5qfql + 5qlql + Bqfq 3 + 5q 2 2 q 3 - ql - 5q 3 q 2 4 

+ 5q 3 qj + hq 3 q\ + 5q 2 ql + 5q 2 2 q 3 + 5g 3 2 g 4 3 - qfj (32) 

b) qi + q 2 > q ?J + g 4 and g x + g 4 < q 2 + g 3 , 



£>i = Qa (33) 
1 

3' 



^2 = \q\ , (34) 



30 

c) q\ + q 2 < q-i + <? 4 and q 1 + q A < q 2 + q 3 , 



D3 = fa? + 5g 2 2 + 5g 3 2 - qj) (35) 



A = 2 (<?i + 92 + ?4 - 93) (36) 



7 In the general case there are two other possibilities 91 > 92 + 93 + 94 and q% > 91+92+94- 
However they are unphysical and give zero answer for all integrals, D\ = D 2 = D 3 = 0. 
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Do 



1 

12' 



-(gi + g 2 ) 3 - 2gf + 2g 4 3 + 3(g a + g 2 )(g 3 2 + <&)) 



(37) 



-D3 is the same as eq. fl32|) but with the change of variables 1 <-> 3 and 2^4. 



d) gi + g 2 < ?3 + g 4 and gi + g 4 > g 2 + g 3 , 



D 2 



6 



£>i = g 2 

g 2 (3g 3 + 3qj - 3g? - qfj 



Ds = ^q! (5g? + 5g 3 2 + 5g 



ql 



(38) 
(39) 

(40) 



Now, the energy 5-function can be trivially integrated away and we are left with 
the following two-dimensional integral upon energies of incoming particles: 



df(E u t) 



1 



PW1PW2 



E' 2 



(41) 



dt 64 7 r 3 £ 1 p 1 

where E 2 = E[ + E 2 — E\ , p 2 = J E\ — m 2 and the summation in eq. (O) is made over 



all Dk contributing to the squared matrix elements in eqs. (|22| , |23[ ) . In the integration 
one should take into account the ^-functions corresponding to the )-d) as well 

as the ^-function related to the energy conservation, 9(Es + £4 — E\ — m 2 ). 

We can further check our calculations of the D functions, using the following trick. 



In the case m\ + m\ = m\ + m\ the l.h.s. of the eq. (|22| ) can be rewritten as: 



Ki(qi^)(qsX) = K^q,^) 2 . (42) 
We can now integrate the r.h.s. of eq. (|42"D over angles. Consequently we obtain 



integrals in the form of eqs. fl25| , p6|) and one more integral: 



£>4(3,4) 



4g|g| f°° d\ 



A 2 



7T 



. . , . . s [ . . . 2 / sin(Ag 3 )\ 

sm(Agi) sm(Ag 2 ) sm(Ag 3 ) - — cos(Ag 3 ) 

Ag 3 \ Ag 3 J 



sm(Ag 4 ) - — cos(Ag 4 ) 

Ag 4 V Ag 4 J 



(43) 
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The equivalence of the two approaches requires that the following two conditions 
should be fulfilled: 

D 2 (l, 2) ee ti + ti-ti-ti Dl + D 2 ( 3 , 4) , (44) 



and 



D 3 ee ^±^1— ^£ 2 (3, 4) + £> 4 ( 3 , 4) . (45) 



We have checked that this is indeed true, and thus assured that the D functions have 
been calculated correctly. 



B The collision integrals in the Boltzmann approx- 
imation 

We have also checked the correctness of our calculations in the following way. In 
the Boltzmann approximation we have calculated everything twice, first letting the 
computer solve the equations directly, and second, by using different exact expressions 
in a half of the collision integral. In the Boltzmann approximation we reduced the 
first half of the collision integrals to the following one dimensional integrals over 
momentum. 

For the reactions v e v e — > v e v e and z/ e 77 e — > v e V e the first half of the collision integral 

is: 



I- Coll 



2 5 G 2 F 



d 3 q 2 d 3 q 3 d 3 q 4 



(27r) 4 5 4 ( gi + q 2 - q, - 94) • F(f) ■ 



where F(f) = 
ing particles. 



2E 1 J (2ir) 9 8E 2 E 3 E 4 
■ 2(g 1 -g 2 ) 2 + 4(g 1 -g 4 ) 2 

— f(in i)f(in 2) with f in being the distribution functions for the incom- 
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For the first 5 reactions in Table p] (including neutrinos only) we obtain the result: 

2 



j first 5 
l coll 



%G 2 F 

3^ 



dp2 P 2 



-f Ue (l) ^ e (2) + -4 i (2)) 



(46) 



The reaction v e v e 



e 1 e gives: 



TVeVe 

1 Coll 



TV 



Et / d£dp 2 pl\ l 



Am 2 



(l-0M-/,e(l)^(2)] 



i 2 2 

-tf L + g fo(l--±)+2g L g R -e- 

o S S 



(47) 



where s = 2EiE 2 (l—^) is expressed trough the angle £ = cos 8i 2 between the incoming 
particles. 



The reactions — > z/ e e ± give: 



coii 



G 2 

— f / G?£dp 2 P2 



7T 



1 - 



E\E 2 



(P 1 -P 2 ) 2 [-U e (l)fe(2)} 



1 -(gl + g 2 R )(A + ^(l + ^))-2g L g R ^ 
3 s s s 



(48) 



where s = 2EiE 2 {l — /3£)+m 2 and is the velocity of the incoming electron (positron). 
As always we have defined g^ = 1/2 + sin 2 6 W and g R = sin 2 6 W . The angle integrals 
are easily done, giving for the v e e^ case: 

P2 



1 coll 



167T 2 Jo 



dp2 T^F 

-C/l-C/2 



[-f Ue (l)f e (2)]m 



1 10m 11m 2 4m 3 



2 u 2 
+AgL9R I 3^ ~ 3u + ~2 + 3 log ^ 



6 + 9 



log(u) 



(49) 



where it = s/m 2 . 

Similarly for the interactions of muonic neutrinos with themselves or other neu- 
trinos, which are the first 5 reactions in Table ^ we get: 



j first 5 Vp 
'-coll 



%G 2 F 



Ex Jdp 2 p\ [-/„„(!) (4/„„(2) + /.j2) 



(50) 



9vr 2 

For the interactions of with electrons we obtain the same expression as for u, 
except for the change g^ — > gi = g^ — 1 and the substitution like f Ue — > /„ . 
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Process 


2~ 5 G F 2 S\A\ 2 


v e + v e 


— ► v e + z/ e 


4(»i • P4)(P2 • P3) 


Ve + V e 


— *■ + 


2(pi -p 2 )(P3 -P4) 


V„ 4- Z/„ 


— > V.i \ -\- V..I \ 


(pi -P 4 )(P2 'Pa) 


+ ^V(r) 




(pi -p 4 )(p 2 -P 3 ) 


+ ^(r) 


— »• V e + l^( r ) 


(Pl -P 2 )(P3 -P4) 


Z/ e + V e 


-> e+ + e~ 


4[M(pi -P4)(P2 - Pa) 

+0fl(Pl -Ps)(P2 -P4) 

+9LgRm 2 e (p 1 -p 2 )] 










v e + e~ 


-> z/ e + e" 


4[^£(P1 -P 2 )(P3 -P4) 

-P4)(P2 - Pa) 

-gLgRmKpi -p 3 )] 










Ve + e + 


-> z/ e + e + 


4[#|(Pl -P 2 )(P3 -P4) 

+#£(Pi -P4)(P2 -Ps) 
-9L9Rml(pi -p 3 )] 











Table 1: The matrix elements for various processes with electronic neutrinos; here 
<?l = 1/2 + sin 2 ^ and g R = sin 2 



Process 


2- 5 G F 2 5|A| 2 


^ + ^ M 




4(pi -p 4 )(p2 -Ps) 


^ + ^ M 


-»• ^ + ^ 


2(P1 -P 2 )(P3 -P4) 


^ + ^ M 


— > Z/ e ( T ) + P e ( T ) 


(pi -p 4 )(p 2 -P 3 ) 


+ V e (r) 


— >• ^ + ^e(r) 


(pi -p 4 )(p 2 -P 3 ) 


+ V e (r) 


— >• ^ + ^e(r) 


(Pl -P 2 )(P3 -P4) 


V(l + Vfi 


-»• e+ + e~ 


4[(<?£(pi -P4)(p 2 - Ps) 

+g 2 R (pi -pz){p2-pa) 

+9L9Rm 2 e (pi ■ p 2 )] 










V? + e" 


-»• ^ + e~ 


4[<?£(pi -P 2 )(P3 -Pi) 
+9r(Pi -P4)(P2 - Ps) 
-~9l9r^I{P\ - Ps)] 










^ + e + 


-> ^ + e+ 


4[#|(P1 -P 2 )(P3 -P4) 
+ <?£bl -P4)(P2 "Ps) 

-9L9Rm 2 e(Pi -Ps)] 











Table 2: The matrix elements for various processes with muonic or tau neutrinos; 
here ql = 9l — 1 = — 1/2 + sin 2 9w and g# = sin 2 6w- 
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Program 


T y /T v 






Full distribution f(x,y) 
momentum points 100 
initial time = 0.2 
Bulirsch-Stoer method U 


1.3996 


0.80 % 


0.26 % 


Full distribution f(x,y) 
momentum points 100 
initial time Xi n = 0.1 
Simple time evolution 


1.3996 


0.79 % 


0.25 % 


Full distribution f(x,y) 
momentum points 200 
initial time Xi n = 0.1 
Simple time evolution 


1.3994 


0.83 % 


0.33 % 


Correction S(x, y) 
momentum points 200 
initial time X{ n = 0.1 
Simple time evolution 


1.3991 


0.94 % 


0.40 % 


Same, but assuming 
entropy conservation 
(photon cooling due to 
i/e-interaction is neglected ) 


1.4010 


1.13 % 


0.53 % 



Table 3: Comparison of the results of different ways of calculation. All the results 
are for FD case. 
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Figure Captions: 



Fig. 1 The ratio T^/T v as a function of the dimensionless "time" x = ma for the 
cases of Fermi-Dirac and Maxwell-Boltzmann statistics. The value x — 1 corresponds 
to T v = lMeV. 

Fig. 2 The relative distortion of the total neutrino energy density Sp/p eq (for the 
definition see eq. (|TTD) as a function of x = ma for the case of "entropy conservation" 
(upper curve) and energy conservation according to eq. (|J) (lower curve). In the first 
case we neglect the energy exchange between neutrinos and e which is equivalent to 
the requirement of entropy conservation in the electromagnetic plasma. 

Fig. 3 The "time" evolution of the correction to the neutrino distribution functions 
8j — (fuj ~ ft q )l fu 9 f° r several fixed values of momentum y = 3,5, 7. In Fig. 3a and 
3b we present respectively the results for electronic neutrinos, and for muonic (tau) 
ones. 

Fig. 4 The evolution of the non-equilibrium corrections to the distribution 
functions for the momentum value y = 5 for the electronic (solid curve) and muonic 
(tau) (dotted curve) neutrinos in the cases of FD and MB statistics. 

Fig. 5 The distortion of the neutrino spectra 8j = (f v . — fl q )/ fl q as functions of 
the dimensionless momentum y at the final "time" x = 60. The dashed lines a and c 
correspond to Maxwell-Boltzmann statistics, while the solid lines b and d correspond 
to Fermi-Dirac statistics. The upper curves a and b are for electronic neutrinos, while 
the lower curves c and d are for muonic (tau) neutrinos. All the curves can be well 
approximated by a second order polynomial in y like 5 = Ay(y — B). 
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